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ABSTRACT 

We describe an updated version of our hybrid A^-body-coagulation code for planet 
formation. In addition to the features of our 2006-2008 code, our treatment now includes 
algorithms for the ID evolution of the viscous disk, the accretion of small particles in 
planetary atmospheres, gas accretion onto massive cores, and the response of A^-bodies 
to the gravitational potential of the gaseous disk and the swarm of planetesimals. To 
validate the N-body portion of the algorithm, we use a battery of tests in planetary 
dynamics. As a first application of the complete code, we consider the evolution of 
Pluto-mass planetesimals in a swarm of 0.1-1 cm pebbles. In a typical evolution time of 
1-3 Myr, our calculations transform 0.01-0.1 disks of gas and dust into planetary 
systems containing super-Earths, Saturns, and Jupiters. Low mass planets form more 
often than massive planets; disks with smaller a form more massive planets than disks 
with larger a. For Jupiter-mass planets, masses of solid cores are 10-100 Mq. 

Subject headings: planetary systems - solar system: formation - stars: formation - 
circumstellar matter 



1. INTRODUCTION 



Gas giant planets form in gaseous disks surrounding young stars. In the 'core accretion' 



rapidly accrete gas (e.g.. 


Poll 


ack et al. 


1996 


Alibert et al. 


2005; 


Ida k. Lin 


2005; 


Chabrier et al. 


2007 




Lissauer & Stevenson 


2007: 


Dodson-Robinson et al. 


200e 


; 


D'Aneelo et al. 201ol). As the 



planets grow, viscous transport and photoevapo ration also remove material from the disk (e.g.. 



Alexander fc Armitage 



2009 



Owen et al. I l2010l ). Eventually, these processes exhaust the disk. 



leaving behind a young planetary system. 



Observations place many constraints on this theory. Although nearly all of the youngest stars 
are surrounded by massive disks of gas and dust, very few stars with ages of 3-10 Myr have dusty 
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disks (e.g., iHaisch. Lada. Ladall200ll : ICurrie et al. II2009I : iKennedv fc Kenvon 1 12003 : iMamaiek 
2009 ) . Thus, gas giant planets must form in 1-3 Myr. Stars with ages of ^ 1 Myr have typi 



cal disk masses o f 0.001-0.1 Mq and disk radii of 10-200 AU (e.g., [Andrews Wilhams I boOSl : 
Isella et al. II2009I ). setting the local environment where planets form. For stars with ages of ~ 1 
Gyr, the frequency of ice/ gas giant planets around low mass stars is ~ 20%-40% (jCumming et al 



2008 



Gould et al. Il2010l ). Although many known exoplanets are gas giants with masses compa- 



rable to Jupiter, lower mass planets ou t number Jupiters by factors of ^ 10 (IMayor et al. 



Borucki et al. 



20101 : 



Gould et al. 



2010 



Holman et al. 



201C 



Howard et al. 



200S 



2ninl v 



Despite the complexity of the core accretion theory, clusters of computers are now capable of 
completing an end-to-end simulation of planet formation in an evolving gaseous disk. Our goal is to 
build this simulation and to test the core accretion theory. Over the pas t decade, we have developed 



a hybrid, multiannulus A^- b ody-coagulation code for p lanet formation (iKenyon fc Bromley 



Bromlev fc Kenvon I l2006l : iKenvon Sz Bromlev I l2008l ) . The centerpieces of our approach are a 



2n04a 



multiannulus coagulation code - which treats the growth and dynamical evolution of small objects 
in 2D using a set of statistical algorithms - and an A^-body code - which follows the 3D trajectories 
of massive objects orbiting a central star. With additional software to combine the two codes into 
a unified whole, we have considered the forma.tion of terrestrial planets (IKenvon &: Bromlev 1120061 ) 



and Kuiper belt objects (IKenvon &: Bromlev 



2004c 



Kenvon et al. 200S ) around the Sun and the 



formation and evolution of debris disks around other stars (jKenvon fc: Bromlev Il2004al . 2004b, 2008, 
2010). 



Treating the formation of gas giant planets with our 2008 code (jKenyon &: Bromley I l2008l ) 
requires algorithms for additional physical processes. As described in §2, we added (i) a ID solu- 



tion f or the evolution of a viscous disk on a grid extending from 0.2 AU to 10'* AU ([Chambers 
2003 : 1 Alexander &: Arrnitage II2009I). fii) algo rithms for accretion of small particles encountering a 



planetary atmosphere (jinaba &: Ikoma I l2003l ) and for accretion of gas onto massive cores, (iii) a 
prescription for the evolution of the radius and luminosity of a gas giant planet accreting material 
from a disk, and (iv) new code to treat the gravitational potential of the gaseous disk and the swarm 
of planetesimals in the A^-body code. Tests of these additions demonstrate th at our complete code 
reproduces the si mulations of visc o us di sks in I Alexander &: Armitaee I (1200911 . accretion rates for 



small particles in llnaba &: 



of A^-body code tests from 



koma I (I2OO3II and fo r 



gas m 



D'Angelo et al. I (|20ld ). and the battery 



Bromley &: Kenyon I (120061 ). In ^3 , we show that our cod e also repro 



duce s results for planets migrating through planetesimal disks (jHahn &: Malhotralll999l : iKirsh et al 
200d ). 



As a first application to gas giant planet formation, we consider whether ensembles of Pluto- 
mass cores embedded in gaseous disks with masses of 0.01-0.1 Mq can become gas giant planets. 
Our results in §4 demonstrate that a disk of Pluto-mass objects will not form gas giants in 10 Myr. 
However, a few Pluto-mass seeds embedded in a disk of 0.1-1 cm pebbles can evolve into a planetary 
system with super-Earths, Saturns, and Jupiters. If we neglect orbital migration, ice and gas giants 
form at semimajor axes of 1-100 AU on timescales of 1-3 Myr. Thus, these calculations match some 
of the observed properties of exoplanets without violating the constraints imposed by observations 
of gaseous disks around the youngest stars. 
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THE MODEL 



2.1. Disk evolution 



Disk evolution sets the context for our planet formation calculations. Viscous processes within 
the disk transport mass inwards onto the central star and angular momentum outwards into the 
surrounding molecular cloud. Heating from the central star modifies the temperature structure 
within the disk and removes material from the upper layers of the disk. All of these processes 
change disk structure on timescales comparable to the growth time for planetesimals. Thus, planets 
form in a rapidly evolving disk. 

For a disk with surface density S and viscosity conservation of ang ular momentum and en- 



ergy l eads to a non-l inear diffusion equation for the time evolution of S (e.g.. lLvnden-Bell &: Pringle 
I974I : |Pringle1ll98lh . 



as 

'dt 



-1 



d_ 
'dR 



/9S 
Vdt 



(1) 



ext 



where R is the radial distance from the central star and t is the time. The first term is the 
change in S from viscous evolu tion; the second term i s the change in S from other processes 



such as photoevaporation (e.g.. [Alexander et al. I l2006l : iQwen et al. I I2OIOI ) or planet formation 



(e.g., I Alexander &: Armitage I l2009l ) . The viscosity is = aCgH, where c<j is the sound speed. 



H is the vertical scale height of the disk, and a is the viscosity parameter. The sound speed is 
= '^kTd/ iiViiH, where 7 is the ratio of specific heats, k is Boltzman's constant, Td is the midplane 
temperature of the disk, /i is the mean molecular weight, and mu is the mass of a hydrogen atom. 
The scale height of the disk is H = c^il^^, where i7 = GM^,/ is the angular velocity. 

There are two approaches to solving equation ([T|). Analytic solutions adopt a constant mass 
flow rate M = Svrz^S through the disk, approximate a vertical structure, and solve directly for 
T,{R,t), T(i{R,t), and other physical variables. Numerical solutions either adopt or solve for the 
vertical structure and then solve for the time va riation of M and other physical variables (e.g.. 



Hueso fc Guillot 11200,4 Mitchell fc Stewart ll201fll V 



Here, we con sider planet formation using both approaches. For the analytic solution, we follow 
Chambers I (j2009l ). who derives an elegant time dependent model for a viscous disk irradiated by 
a central star. For the numerical solution, we assume that the midplane temperature is derived 
from the sum of the energy generated by viscous dissipation (subscript "V") and the energy from 
irradiation (subscripit "I") 

(2) 



The viscous temperature is 



rp4 



(3) 



wher e n is the opacity and a is the Stefan-Boltzmann constant (e.g.. lRuden &: Lin lll986l : lRuden k. Pollack 



I991I ). With 1/ = aCg Q.-'^ and t2 = (27a/64cr) {'jk/iimH) kOS^, the viscous temperature is 



Ty = t2Td. 



(4) 
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If the disk is vertically isothermal, the irradiation temperature is Tf{R) = (0/2)(i?^/i?)^T^, 
where R^, and are the radius and effective temperature of the central star and 



2,-K \ R 



d{H/R) 
dR 



([Chiang &: Goldreich Ill997l ). Thus, the irradiation temperature is 



2 

3^ 



R. 
R 



+ 



H /A 



2R\ R 



dlnH 
d In R 



(5) 



(6) 



Following IChiang k Goldreich I (|l997l l and lHueso fc Guillot I (|2005l l. we set d In H/d In R = 9/7. 
With H = c,0-\ we set to = {2T^ /^t^){R^/ Rf and ti = {R^/Rf [IRQ)-^ (jk/^imH)^^^ T^- The 
irradiation temperature is then 



Tf = to + hT^ 



1/2 



(7) 



Altho ugh viscous disks are not vertically isothermal (|Ruden &: Pollack I Il99ll : iD'Alessio et al. 



19981 ). this approach yields a reasonable approximation to the actual disk structure. 



Because Ty and T/ are functions of the midplane temperature, we solve equation ([2]) with a 
Newton- Raphson technique. Using equations ([4]) and ([7]), we re- write equation ([2]) as 



f{Td 



Adopting an initial ~ ^2^'^ °^ '^d ~ ^i' ' ; the derivative 



T| - (to + tiT]/' + t^Ta) 

2/7 







4TS - —T 



tl -1/2 



t2 



allows us to compute 



5Td = f 

and yields a converged to a part in 10® in 2-3 iterations. 



(diy 



(8) 



(9) 



(10) 



In the inner disk, the temp erature is of t en ho t enough to vaporize dust grains. To account for 
the change in opacity, we follow [Chambers I (j2009l ) and assume an opacity 



(11) 



1991 



SteDinskilll998l l. For 



with n = -14 in regions with Ta > = 1380 K (jRuden fc Pollack 
simplicity, we assume k = kq when < T^. 

To solve for the time evolution of S, we use an explicit techniq ue with N annuli on a grid 

, 1982). We adopt an initial 

(12) 



extending from Xin to Xout where x = 2 R^/"^ ( Bath &: Pringle 
surface density 

So 



1981 



2'kRRq 

where M^ n is the initial dis k mass and Rq is the initial disk radius (jHartmann et al 



1998; 



Alexander &: Armitage I l2009l ) . For the thermodynamic variables, we adopt 7 = 7/5 and ^ 
2.4. 
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Figure [T] compares analytic and numerical results for a disk with a = 10~^, M^^ = 0.04 M©, 
and Rq = 10 AU surrounding a star with = 1 Mq. The numerical solution tracks the analytic 
model well. At early times, the surface density declines steeply in the inner disk (where dust grains 
evaporate) and more slowly in the outer disk (where viscous transport dominates). At late times, 
irradiation dominates the energy budget; the surface density then falls more steeply with radius. 

Figure [2] compares the evolution of the disk mass and accretion rate at the inner edge of the 
disk. In both solutions, the disk mass declines by a factor of roughly two in 0.1 Myr, a factor of 
roughly four in 1 Myr, and a factor of roughly ten in 10 Myr. Over the same period, the mass 
accretion rate onto the central star declines by roughly four orders of magnitude. 

In addition to visc o us ev olution, we consider photoevaporation of disk material. Following 



Alexander Sz Armitage I (|2007l ). we assume high energy photons from the central star ionize the 
disk and drive a wind. As long as the inner disk is optically thick, recombination powers a 'diffuse 
wind,' where the change in surface density is 

.as. 



g^-)diffuse = -2no{R)uo{R)nmH ■ (13) 



In this expression. 



and 



/ p N 0.2457 

uoiR) = 0.3423 c,,^ (^-^ - 0.1 j e-^-^^i^ (r/r.-o.i) ^-^^^ 



The recombination coefficient is = 2.58 x 10 cm^ s ^ ( Cox Il2000l ). The gravitational radius 
is Rg = GMi,/c1^, where Cs^w = 10 km s~^ is the sound speed in t he wind. For typica l stars, the 



luminosity of H- ionizing photons, is ~ 

1040_io42 g-i. However. lOwen et al. I (120101 ;) show that 
X-rays drive a more powerful wind than lower energy photons. Here, we assume that the mass loss 
rate in the wind, M^i^d = / 2ttR dTi/dt dR, is a free parameter and consider disk evolution for a 
range in M^ind- This approach is equivalent to assuming a range in <I>^. 

As the disk evolves, the diffuse wind removes more and more material throughout the disk. 
Eventually, the inner disk becomes optically thin to ionizing photons. Ionization then drives a 
'direct wind.' The change in surface density with time is then much larger, 

(-),_. = 0.47 M mu c,. (4,„^(^^)^3^(,) J (^J > ^ > ^^^^ (16) 

where Rin{t) is the inner edge of the optically thick portion of the disk. Inside Rin{t), the mass 
loss rate from the disk is zero. 

For mass loss rates M^i^d ~ lO^"*^" — 10~^ Mq yr~-^, the surface density evolution of a viscous 
disk with photoevaporation follows the evolution in Figure [1] for ~ 10^ — 10^ yr. As the accretion 
rate through the disk drops, mass loss from the wind becomes more and more important. Once the 
inner disk becomes optically thin, the direct wind rapidly empties the inner disk of material. As 
the system evolves, the size of this inner 'hole' grows from Rh ^ Rg 3 AU to Rh ~ 30-100 AU 
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in 0.01-0.1 Myr (see also [Alexander &: Armitage II2009I : lOwen et al. II2OIOI ). A few thousand years 
later, the gaseous disk is gone. 

Figure [3] shows the variation of disk mass and mass accretion rate for the disk in Figure [2] 
and several mass loss rates. In general, the wind starts to empty the disk when the mass accretiori 
rate through the d isk falls below the mass loss rate in the wind (see [Alexander &: Armitage II2009I : 
Owen et al. II2OIOI . and references therein). For very low mas s loss rates, M^i,inA = 10~^^ M fr, yr~^, 



the wind evaporates the disk on timescales of 10 Myr (see also [Alexander &: Armitage 1120071 . 2009). 
As the mass loss rate from photoevaporation grows, the disk evaporates on shorter and shorter 
timescales. For Mwind = 10~^ Mq yr~^, the disk disappears on timescales shorter than 1 Myr. 

These evolutionary models capture the main observable features of disks around pre-main 
sequence stars with ages of 1-10 Myr. For 1 Myr-old young stars, disk masses of ^ 0. 001-0.1 
Mf7^ and mass accretio n rates of 1 — 100 x lO"^*^ M f^ yr~^ are common ( Gullbring et al 



Hartmann et al. 



(Hartmann et al. 



19981 : 



19981: [Andrews Sz Wilhams I l2005l. 2 0071. Typical disk lifetimes are 1-3 Myr 



19981 : iHaisch. Lada, &: Ladall200lh. Few pre-main sequence stars h ave ages 



larger than ~ 10 Myr (jCurrie et al. I l2007l : iMamaiek I l2009l : iKennedv k Kenvon 1 12009| ). Thus, 



these calculations yield reasonable physical conditions for planet formation. 

In our previous calculations, we set the surface density of the gaseous disk as S = So2;ma~"e~*/*s , 
with scaling factor, n as a constar it power law, and tg as a constant gas depletion time. In 

our new approach, adopting an analytic ([Chambers II2009I ) or a numerical (eq. [T|) solution to the 
radial diffusion equation has several important advantages for little computational effort. 



• Following the evolution of the mass accretion rate onto the central star enables robust com- 
parisons between observations of accretion in pre-main sequence stars with the timing of 
planet formation throughout the disk. 

• Treating the evolution of photoevaporation allow s us to make links between t he so-called 
transition disks and the formation of planets (e.g.. [Alexander &: Armitage II2009I ). 



Track ing the time evolution of the radial position of the snow line (e.g., [Kennedy Kenvon 



20081 ) in a physically self-consistent fashion gives us a way to include changes in the compo- 



sition of planetesimals with time. 

Calculating the radial expansion of t he disk with time provides a way to compare the observed 
properties of protostellar disks (e.g.. [Andrews &: Williams II2005I . 2007, 2009) with the derived 
properties of planet-forming disks. 



It is straightforward to derive the properties of power-law disks that yield results similar to 
those of our new calculations. In analytic and numerical solutions to the diffusio i i equa tion, the 
viscous time ti, ~ E? /Zv is roughly the gas depletion time tg. For the [Chambers I (j2009l ) analytic 
approach, ty ^ A Myr (a/10~^); for our numerical solution, tj^ ~ 1 Myr (a/10~^). In previous 
calculations, w e adopt tg = 10 Myr; t hus, n ew calculations with a = 1 — 4x 10~^ roughly correspond 



to the disks in lKenvon Bromlev I (|2008l . 2009, 2010). 
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The new calculations generally yield shallower power law slopes - n = 0.6 for [Chambers 

(|2009l ) and n = 1 for our numerical solution - than the n = 1.0-1. 5 assumed in our previous 
studies. Because the timescale for planet formation is t oc a"^^'^ (see iKenvon Bromley I l2008l . 
and references therein), planets form relatively faster in the outer disk in these new calculations. 
However, the range of dis k masses considered here, M^n = 0.01-0.1 Mq, overlaps the M^q = 
0.0003-0.25 Mq adopted in lKenyon Sz Bromley I (|2010l ). Thus, calculations with similar initial disk 
masses and similar initial size distributions of planetesimals will yield similar timescales for the 
formation of the first oligarchs. 



In 



As a concrete example, w e compare several disk surface density distributions quantitatively. 



Kenvon Sz Bromlev I (j2008l . 2009, 2010), we often adopted an initial surface density distribution 
of solid material, S = 30 g cm~^ Xm (a / 1 AU)~^/'^. At 5 AU, our numerical solutions for disks 
with Md = 0.1 Mq and Rdo = 30 AU have the s ame surface densit y as power law disks with Xm 



= 3. For the same initial disk mass and radius, the [Chambers I (j2009l ) analytic model has the same 
surface density as power law disks with Xm ~ 2.5. For identical initial conditions and equivalent 
viscous timescales, these three surface density distributions yield similar evolutionary times for the 
growth of protoplanets. 



2.2. Coagulation Code 



Kenvon fc Bromlev I (j2008l ) describe the main details of the coagulation code we use to calculate 
the growth of small solid particles (planetesimals) into planets. Briefly, we perform calculations on a 
cylindrical grid with inner radius Rin and outer radius Rout- The model grid contains concentric 
annuli with widths 5Ri centered at radii Ri. Calculations begin with a mass distribution n{mik) 
of planetesimals with horizontal and vertical velocities hik{t) and Vik{t) relative to a circular orbit. 
The horizontal velocity is related to the orbital eccentricity, e^^(t) = 1.6 {hik{t) /VK,if' , where Vk^i 
is the circular orbital velocity in annulus i. The orbital inclination depends on the vertical velocity. 



ik 



it) 



sm 



K,i) 



The mass and velocity distributions of planetesimals evolve in time due to inelastic collisions, 
drag forces, and gravitational forces. The collision rate is navfg, where n is the number density 
of objects, a is the geometric cross-section, v is the relative velocity , and fg is the gravitational 



focusing factor (jSafronov 



Weidenschilling et al. 119971: 



1969 



Lissauer 



1987: Spaute et al. 



1991 



Wetherill Stewart 



uu lll998l : lK^ivov et al. Ibood : iThebault fc Augereau 



1993 



Lohne et al. 



Kenvon &: 



2007 



20081 : Kenvon Sz Bromlev 20081 ). The collision outcome depends on the ratio of the 



collision energy needed to eject half the mass of a pair of colliding planetesimals to the center 
of mass collision energy Qc- If mi and m2 are the masses of two colliding planetesimals, the mass 
of the merged planetesimal is 

m = mi + m2 — mgj , (17) 



where the mass of debris ejected in a collision is 



m. 



0.5 (mi + m2) 



Qc 
Q*, 



9/8 



(18) 
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This approach ahows us to derive ejected masses for catastrophic col l isions with Qr ^ Q*[ and for 



Tanaka et al. 



1996 



Kobavashi &: Tanaka 



cratering collis i ons w i th Qr <^ (see also metherill &: Stewar t 1993: IWilliams Sz Wetherill 



Stern k Colwell I Il997l : iKenvon Luu I il999i : lO'Brien k Greenberg 



1994 



20031 : 



2010I'). Consistent with numerical simulations of c ollision outcomes (e.g. 



(19) 



Benz Sz Asphaug Ill999l : iLeinhardt et al. 1120081 : iLeinhardt fc Stewart II2009I ). we set 

Q* = Q,r^^ + Qgp^^ 

where Qb^^^ is the bulk component of the binding energy, QgPgr^^ is the gravity component of the 
binding energy, r is the radius of a planetesimal, and pg is the mass density of a planetesimal. 



To compute the evolution of the velocity distribution, we include collisional damping from 
inelastic collisions and gravitational intera ctions. For ine lastic and elastic collisions, we adop t 
the statistical, Fokker-Planck approaches of lOhtsuki I (jl992l ) and lOhtsuki. Stewart, k Ida I (|2002l ). 
which treat pairwise in teractions (e.g., dynam ical frict ion and viscous stirri n g) be tween all objects 

we add terms 
see also 



(2001 



in all annuli (see also IStewart k Ida I l2000l ) . As in iKenyon k Bromley 
to treat the probability th at objects in annulus i interact with objects in annulus j 
Kenvon k Bromlev l[2004bl . 2008). 



2.2.1. Small particles 



In most of our previous calculations, we calculate the evolution of particles with radii larger 



than the 'stop ping radius. 



0.5-2 m at 5-10 AU (jAdachi et al. Ill976l : IWeidenschilhng 



1977 



Rafikov I l2004l ). Although they are subject to gas drag, these particles are not well-coupled to 
the gas. Here, we also co nsider the evolution of smaller particl es entrained by the gas (see also 
Kenyon k Bromley II2009I ). Following lYoudin k Lithwick I ([20071), we derive the Stokes number 

rpgQ 



St 



PCs 



where p = T,/2H is the gas density. The vertical scale height of small particles is then 



H 



St<a 
St > a 



(20) 



(21) 



We assume that entrained small particles have vertical velocity v 
h = 1.6v. 



HgQ and horizontal velocity 



Estimating accretion rates of small particles by embedded protoplanets is a challenging prob- 
lem. For particles with St ^ 1, accretion rates depend on the strength of the p lanet's gravity 
relative to forces that couple particles to the gas. Recently, lOrmel k Klahr I (j20ld ) began to ad- 
dress this issue with an analysis of interactions between protoplanets and particles (loosely) coupled 
to the gas. After identifying three classes of encounters, they derive growth rates as a function of 
m, St, and the local properties of the gas. For massive protoplanets with r > 1000 km, they derive 
short growth times, ~ 10^ yr, for particles with 10^^ < St < 10. 

Here, we adopt a simple prescription for accretion of small particles with St < 1. Small 
protoplanets with r < 100 km do not have strong enough gravitational fields to wrest small particles 
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from the gas. Thus, w e assume small protoplanets cannot accrete particles with St < 1. The 
Ormel &: Klahr I ( 20101 ) results suggest large protoplanets accrete particles with St < 10^^ on long 
timescales, > 1 Myr. Thus, we also assume protoplanets of any size cannot accrete these particles. 
For particles with St > 10~^, we assume accretion proceeds as in our standard formalism with 
collision rates navfg, where the cross-section a includes the enhanced radius of the protoplanet 
di scussed in ^2.2.2. Our a pproach generally yields growth rates a factor of 1-10 smaller than those 
of lOrmel &: Klahr I (j2010l ). Thus, our formation times are longer than timescales derived from a 
more rigorous approach. For the calculations in this paper, however, this accuracy is sufficient 
to explore the impact of disk mass and viscosity on the formation of gas giant planets. In future 
papers, a more comprehensive treatment of small particle accretion and gas drag will yield a better 
understanding of the role of the local properties of the gas on planet formation. 

Aside from enabling protoplanets to accrete smaller particles at large rates, including small 
particles with r < 1 m allows us to derive a more accurate calculation for the t i me ev olution 
of debris from the collisional cascade. Compared to results in iKenyon &: Bromley I (l2008l . 2010), 
estimates for the abundance of very small grains with r ~ 1-10 fiui now rely on extrapolation over 3 
orders of magnitude in radius instead of 6. Thus, this addition yields m ore accurate predictions for 



the variation of grain emission as a function of time (see, for example, iKenvon &: Bromlev 
2010). 



200S 



2.2.2. Planetary atmospheres 



As planets grow, they acquire gaseous atmospheres. Initially, the radius of the atmosphere Ra 
is well- approximated by the smaller of the Bondi radius Rb and the Hill radius Rh- 



Ra 



Rb 



GMr, 



1/3 



Rb < Rh 
Rb > Rh 



(22) 



where Mp is the mass of the planet and a is the semimajor axis of the planet's orbit around the 
central star. For low mass planets, Rh > Rb- Requiring Rb > Rp, where Rp is the radius of 
the planet, an icy planet with a mass density of 1 g cm~^ starts to form an atmosphere when 



Mp > 3 X 10 



25 



Once planets develop a small atmosphere, t hey accrete small particles m ore rapidly (jPodolak et al 



19881 : iKarv et al. Ill993l : 



Inaba &: Ikoma II2003I : iTanigawa &: Ohtsuki Il2010l ). As small particles ap 



proach the planet, atmospheric drag reduces their velocities. S maller velocities i ncrease gravita - 
tional focusing factors, enabling very rapid accretion (see also [Chambers I l2008l : iRafikov I l2010l ). 
Because the extent of the atmosphere depends on the accretion luminosity, calculations need to 
find the right b alance between the exte nt of the atmosphere and the accretion rate (and luminosity) . 
Here, we follow Inaba fc Ikoma (j2003l ) and solve for the structure of an atmosphere in hydrostatic 



equilibrium around each planet with Mp > 3 x 10^^ g. For each smaller planetesimal with mass m, 
we calculate the enhanced radius Re{Mp, m), which replaces the physical radius in our derivation 
of cross-sections and collision rates. 



This approach shortens growth times by factors of 3-10. In 



Chambers 



(j2006l ). protoplanets 
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with atmospheres gro w 3-4 times faster at 5 AU than protoplanets without a tmospheres (see also 
Inaba &: Ikoma To confirm these results, we repeated the calculations of lKenvon &: Bromley 



(j2009l ) for protoplanets with atmospheres. In a power-law disk with Xm = 5, our previous calcu- 
lations yielded 10 protoplanets on timescales of ~ 3 Myr. Repeating these calculations for 
protoplanets with planetary atmospheres results in growth times shorter by factors of 5-10, ~ 0.3- 
0.5 Myr. Several numerical tests suggest that different fragmentation laws produce factor o f 1.5-2 
changes in the grow t h tirn e for protoplanets with atmospheres (see also IChamberd l2006l . 2008). 
Kenyon fc Bromley (j2009l ) derived only 10% to 20% variations in growth times for a similar range 



in fragmentation laws. Thus, growth times are more sensitive to fragmentation when protoplanets 
have atmospheres. 



2.2.3. Gas accretion 



As planets continue to grow, they begin to accrete gas from the disk. Although the minimum 
mass required to accrete gas varies with the accretion rate of the planet, the properties of the disk 



( Mizuno 


1980: 


Stevenson 


1982: 


Ikoma et al. 


2000: 



1-10 Mg 

2010). Once gas accretion begins, the planet mass grows rapidly until it removes most of the gas 
in its vicinity, either by depleting the disk entirely or by opening up a gap between the planet and 
the r est of the disk. The planet then grows mo r e slowly and may start to migrate t hrough the disk 



(e.g., iLin Paoaloizou Ill986l : iLin et al. Ill996l : I Ward Ill997l : IP'Angelo et al. II2OIOI ). 



Here, we approximate the comp licated physics of gas accretion and the evolution of the atmo- 



spher e with a simple formula (see also 
2OO9I ). The growth rates reviewed in 



Veras &i Armitage1l2004 ; Ida &: Liii1l2005 ; Alexander &: Armitage 



D'Angelo et al. I (|201fll ) suggest 



Mo Cs 



-((M-M0)/o-m)2 



(23) 



where Mq is a typical maximum accretion rate, Mq ~ 0.1 Mj (1 Mj is the mass of Jupiter), 
= log Mp, fiQ = log Mq, and am ~ 2/3. This functional form captures the rapid increase in the 
accretion rate from the minimum core mass of 1-10 to larger core masses, the peak accretion 
rate at masses of roughly 10% the mass in Jupiter, and the decline in the accretion rate once 
the planet opens up a gap in the disk. Detailed numerical simulations suggest maximum accretion 



rates of 1 Jupiter mas s every 0.01-0.1 Myr (e.g. iLissauer fc: Stevenson 



2007 



D'Angelo et al. 



201c 



Machida et al. II2OIOI . and references therein). For simplicity, we consider Mq a free parameter; we 
set a rate appropriate for a Jupiter mass gas giant in a disk with S = 100 g cm~^, = 1.78 x 10~^ 
s^^ (5 AU), and Cs = 0.67 km s^^ and use equation (j23p to scale this rate throughout the disk. 

To remove accreted gas from the disk in our numerical simulations of disk evolution, we set 
the change in surface density as 



(-) 



Mr, 



dt ~ 27rRAR ' ^^^^ 
where R is the semimajor axis of the planet and AR is the width of the annulus. Combined with 
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photoevaporation from the central star, this expression yields 



/9S _ j {^)diffuse + {^)gas , Rin ~ R* ,r)r\ 
y r, )ext — \ .qsn \(dT,\ P ^ ' 

'-"^ I ) direct ~r \~gi)gas i ^in 



where Rin is the inner edge of the optically thick portion of the disk. 



Although we use equation (p3|) to derive gas accretion rates for the analytic disk model, we do 
not remove the accreted mass from the disk. To place some limit on the amount of mass accreted 
by gas giants, we halt gas accretion when the total mass in gas giants exceeds the total remaining 
mass in the disk. 



2.2.4- Evolution of R and L of planets 



Throughout this phase, t he radius of the planet Rp depends on the accretion rate and the prop 



erties of the atmosphere f e.g 



2007 



Lissauer fc Stevenson 



Hartmann et al. 



2007 



1997; 



Baraffe et al. 



Papaloizou fc Nelson II2005I : IChabrier et al. 



20091 : iD'Aneelo et al. \\20m . Before the planet 



opens a gap in the disk, accretion is roughly spherical; after the gap forms, the planet accretes 



material from a thin disk- like structure within the planet's Hill sphere (e.g., iNelson et al. Il200d : 
D'Angelo &: Lubow I l2008l ). For planets with masses much larger than Neptune, disk accretion 
provides most of the gas. Thus, to derive a simple estimate for the evolution of Rp, we solve 
the energy equation fo r a polytrope of index n = 1.5 accreting material from a gaseous disk (e.g. 



Hartmann et al. 1119971 ): 



Rr. 



l\ Mr, 



Rp 3J Mp 



3 \GM^/Rp 



(26) 



Here, L is the planet's photospheric luminosity in erg s ^; 1 — r], where r] < 1, measures the fraction 
of the accretion energy radiated away before material hits the planet's photospher^. 



Solving equation (j26p requires a relation for L. When M = 0, planets rese i nble n = 0.5-1 
polyt ropes, with Rp oc and —1/8 < m < 1/10 (e.g., ISaumon et al. 



1996 



Fortnev et al. 



20091'). For this phase, we derive a simple expression from published evolutionary tracks (e.g. 



Burrows et al. 1119971 : ISpiegel et al. 



2010 



Baraffe et al. 112001 2008): 



L ^ 10 Lf^ 



, , (27) 
where Rj is the radius of Jupiter. Accretion energy tends to expand planets considerably (e.g. 



0.25 



Rr) 



15 



Pollack et al. Ill996l : iMarley et al. II2007I : iD'Angelo et al. Il2010l ). To treat this phase, we derive a 
simple expression frora published model at mosphere calculations applied to an n = 1.5 polytrope 



( Saumon et al. 



1996 



Fortnev et al. II2009I . and references therein): 



L ^ 10"'' L, 







Mp 
1 Mj 



R„ 



1 Rj 



(28) 



^To avoid confusion with our use of a for the disk viscosity parameter, we change the a in lHartmann et al. I (1997) 
to rj. 
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When the planet has a radius Rt ~ 1.78(Mp/l Mj) Rj, these relations yield the same luminosity. 



Thus, we use equation ([27|) for R{Mp) > Rt and > and equation ([28]) otherwise. 

Despite its simplicity, this approach yields reasonable results for L{t) and Rp{t). In the example 
of Fig. HHSl the calculation starts with a 10 Mq, 1 Rj planet accreting at a rate of 10~^ Mj yr"^ 
until it reaches a mass of 1 Mj. Before the planet reaches its final mass, the luminosity increases 
exponentially with time. Once accretion stops, the luminosity declines. Peak luminosity and radius 
depend on ry; planets that accrete hotter material (larger rf) expand more than planets that accrete 
colder material (smaller rj). F or 7? ~ 0.3-0.4, our peak luminosities of ^ 10~^ LfD are sim ilar to those 



of more detailed calculations (jPollack et al. Ill996l : iMarlev et al. 1120071 : iD'Angelo et al. Il20ld l. At 



late times, evolution for all 77 converges on a single track for L{t) and Rp{t). This track yields a 
radius of 1.03 Rj at 4.5 Gyr. 

Although this approach is not precise, it serves our purposes. The model yields planetary 
radii with sufficient accuracy (it 10%) to derive the merger rate of growing planets from dynamical 
calculations with our A^-body code. The estimated luminosity is also accura te enough to serve as a 
re asonable starting poin t for more detailed evolutionary calculations, as in iBurrows et al. I (|l997l ) 



or 



Baraffe et al. I (j2003l ). This formalism is also flexible: it is simple to adopt better prescriptions 



for the luminosity or an energy equation with a different polytropic index. 



2.3. A-body Code 



Bromley &: Kenyon I ()2006l ) provide a description of the A-body component of our hybrid code. 
To solve the equations of motion for a set of interacting particles, w e use an adaptive algorithrn with 
sixth-order time ste ps, based either on Richardson extrapolation feromlev fc Kenvon 2006) or a 



symplectic method (lYoshida I Il990l : iKinoshita. Yoshida. &: Nakai 1 ,1991. : .Wisdom &: Holman 



1991 



Saha fc Tremaine Ill992l ). The code calculates gravitational forces by direct summation and evolves 
particles accordingly. In our earlier version, we performed integrations in terms o f phase-space 



variables defined relative to individual Keplerian frames, following the work of Encke (lEncke 



Vasilev Ill982l : IWisdom fc Holmanlll99ll : llda fc Makino Ill992l : iFukushima Ill99fil : [SheferJ 12002); our 



1852 



latest version can track center-of-mass frame coordinates instead. This feature may be useful in 
instances where stellar recoil is important. As before, the code evolves close encounters - including 
mergers - between pairs by solution of Kepler's equations in the pairs' center of mass frame. The 
algorithm also includes changes in eccentricity, inclination, and semimajor axis, derived from the 
Fokker-Planck and gas drag algorithms in the coagulation code. 

The 2006 version and the new version of A-body code can track the force from a dusty or 
gaseous disk. The older code handles only toy models for the gaseous disk, consisting of a power 
law surface density profile and an exponential decay in time. The new code includes the more 
realistic disk potential derived in §2.11 The code represents a disk in annular bins with time- 
varying surface density specified by the physical model. Although the current code does not treat 
the lack of gravitational force from gaps in the disk produced by gas giant planets, this extra force 
is small compared to the forces from the rest of the disk and gas giant planets. Once this feature 
is included, we can then self-consistently track the dynamics of gas giants in an evolving gaseous 
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disk. 

Appendix [XI provides some background and details of our method to treat the gravity of the 
gaseous disk. To summarize, we approximate the disk as axisymmetric with a scale height H, which 
is small compared to the orbital distance {H/a ^ 0.1). The acceleration that the disk produces 
on a planet near the disk midplane is then fast to calculate. When the disk has a power-law 
surface density, a simple analytical expression suffices. If the disk has more complicated radial 
structure ( ^2.ip . then we split it up into constant-density annuli and calculate the contribution of 
each annulus to the acceleration at some point near the disk midplane. We repeat this calculation 
for a set of O(IO^) radial points at the beginning of a simulation, and interpolate with a cubic spline 
thereafter to find the disk acceleration at the location of the objects. As the surface density of 
the disk evolves, the code updates the date for the spline fit. 

In addition to gas, the new version of the A-body code includes t he gravity of non-physical 



object s representing dust or planetesimals. As in the SyMBA code of iDuncan. Levison. fc Lee 



(j 19981 ). we use massless "tracer" particles, which simply respond to the gravitational potential 
produced by the central star and planets. We use the evolution of the tracers to inform the 
coagulation code of the de/dt, di/dt, and da/dt induced by massive planets in the A^-body code. 
We also have a second population of massive "swarm" particles, which have mutual interactions with 
planets but do not interact with either tracers or with each other. The swarm thus consists of super- 
particles which represent the ensemble of small objects that interact with each other statistically 
in the coagulation code and which also interact dynamically with planets in the A-body code. 
Independently of the tracers, we can use the swarm particles to inform the coagulation code of the 
de/dt, di/dt, and da/dt induced by massive planets in the A-body code. More importantly, we 
use the mutual interactions of swarm and A-body particles to calculate the response of planets 
to the surface densi t y and velocity distributions of the swarm. In addition to the tests in §3, 



Bromley fc Kenvon I (j201ll ) describe some applications of the interaction between swarm and A- 



body particles. 

In terms of comparative computational cost, planets are expensive, tracers are cheap, and 
swarm particles are in between. Even without mutual interactions, the swarm particles can be a 
considerable burden in long-term integrations. To lighten this computational load, we evolve the 
star and planets independently of the swarm during a single coarse- grain timestep, allowing the 
code to take as many substeps as necessary. Then, we interpolate the resulting planet trajectories 
to calculate the gravitational forces needed to update the swarm. When we use massive swarm 
particles, we allow the planets to deviate from their interpolated path in response to the smaller 
objects. This approach is valid only when the net forces on the planets from the swarm vary slowly 
in time, and when the our simple third-order interpolation of the planet's trajectory over a coarse 
timestep is realistic. With this approach, there is a risk of lower accuracy for the very rare events 
when some members of the swarm interact with a pair of closely-interacting A^-body particles. 
Our 3rd-order interpolation then does not yield a robust representation of the trajectories of these 
swarm particles. Detailed tests show that reduced accuracy only occurs for the very few swarm 
particles within 1-2 Hill radii of the pair of A-bodies, has a negligible impact on the trajectory 
of the A^-body particles, and does not infiuence the outcome (merger or scattering event) of the 
interaction between the two A-bodies. 
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2.4. A^-body code tests 



We put the new version of the A^-body code through the same diagnostics as in lBromlev &: Kenvon 



pood ). To assess accuracy and dynamic range, we perform long-term orbit integrations of the ma- 
.j or planets, as well as integrations of the planets with t heir masses scaled by a factor of fifty 
(jPuncan. Levison. &: Lee 1 119981 : iBromley &: Kenyon 1 120061 ) . We also track the motion of tight bi- 
naries in orbit around a central star. In one limit, we evolve a close pair of Jupiter-mass objects 
(jPuncan. Levison. Leelll998l ): in another, we simulate a surface-skimming satellite above the 
Earth as it orbits the Sun. We also set up pairs of planets on closely-spaced circular orbits and 
evolve them t o see whether w e can resolve the critical separation that determines whether their 
orbits crosj^ ( Gladman 1993 ). With these tests we verify t hat the new code can replicate the 
results illustrated in Figures 1-4 of iBromlev &: Kenvon I (|2006l ). 



To save computation time, the new code evolves tracers and swarm particles with lower tem- 
poral resolution , than w ith the planets. To verify that the resulting orbits are reasonable, we follow 



Kokubo &: Ida I (|l995l ) in simulating gravitational stirring of a planetesimal swarm by a pair of 
planets. Two 2 x 10^^ g planets are separated by 10 Rh and are embedded in the middle of an 
annulus of 800 2 x 10^'^ g objects, which is 35 Rh wide and centered at a distance of 1 AU from a 
1 Mq star. All objects are initially on circular orbits. We evolve this system treating all particles 
as mutually interacting massive bodies, then with the disk composed of massive swarm particles 
that interact with the planets only, and finally with a massless tracer disk. Figure [6] shows that 



(Kokubo & Ida 


1995; 


Weidenschilline et al. 


1997: 


Kenvon & Bromlev 


2001: 


Bromlev &: Kenvon 


2006 


)• 



The final set of tests for the A^-body code concerns its ability to identify mergers, in particular 
between planets and swarm particles or tracers. We first confirm that minor changes to the code 
m aintain its accuracy in hig h-speed mergers between small massive objects (the "grapefruit test" 



in lBromlev &: Kenvon II2006I ). We then test mergers betweei i planets and tracers or s warrn particles 
explicitly with the more sophisticated method proposed bv iGreenzweig &: Lissaueii ([1990) • We set 
up (i) a 10~^ Mq planet on a circular orbit at 1 AU and (ii) test particles on orbits with eccentricity 
e = 0.007, inclination i = 0.2 degrees, and semimajor axes distributed in two rings (0.977 AU to 
0.991 AU and 1.009 AU to 1.023 AU). We evolve the system through a single close encounter between 
the planet and each test particle to measure the fraction of test particles accreted by the planet. 
With the planet's physical radius s et to 5300 km, we derive a merger fraction of 0.008 ± 0.001 



consistent with previous estimates (jGreenzweig &: Lissauer 
Bromlev &: Kenvon 1120061 ). 



1990 



Duncan. Levison. Lee 



1998 



^Our code yields a value within a few percent of the predicted critical separation, Aocrit ~ 2\^Rh, where Rh is 
the mutual Hill radius (as in eq. [22] but with Mp set to the sum of the planets' masses). 
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MIGRATION 



Interactions between a gaseous or planetesimal disk and planets can lead to signification mi- 
gration, wi th rates of change in sernimajor axis, da/dt, that are important on plaiiet formation 



timescales (iLin &: Papaloizou 



19791: iGoldreich fc Tremaine lllQSd : iLin Papaloizou 



1986 



Ward 



19971 : IVeras &: Armitage I l2004l : iLevison et al. 1120071 ). Migration occurs when a planet perturbs 
a disk gravitationally. The resulting angular momentum exchange between the planet and the 
perturbed di sk causes the planet to move radially inward or outward. Several types of migration 
(jWard 1119971 ) depend on disk and planet properties. When an embedded planet produces a modest 
perturbation in a massive disk (type I migration), the planet responds to torques from the disk 
and migrates inwards. When a planet is massive enough to open a radial gap in the disk (type II 
migration), the planet becomes locked into the global viscous flow and migrates inwards. Type III 
migration (e.g., iPapaloizou et al.l 120071 ) is a fast migration mode associated with smaller planets 



that efficiently exchange angular momentum with co-orbiting material. As with migration in a 
planetesimal disk (see Fig. [T0|) . this mode can produce inward or outward migration. 

Our c ode treats all types of m igration through the gaseous disk using semianalytical approxima- 
tions (e.g.. IPapaloizou et al.ll2007l ). Although we do not include any detailed disk hydrodynamics, 
the code can smoothly change the semimajor axis and other orbital elements of each planet at any 
specfied rate. 

The code includes another migration mechanism, perhaps "Type 0," that involves interaction 
between a massive disk and planets. In an axisymmetric disk that redistributes mass, planets 
migrate as the disk mass changes. For photoevaporation and viscous transport, the disk mass 
and gravitational potential change slowly compared to the local dynamical time. The angular 
momentum of the planet is conserved. By treating the disk and the central star as a point mass 
with mass Me// acting on a planet at semimajor axis a, aM(>ff is conserved. As the disk vanishes, 
the planet migrates outward. We provide a few more details of Type migration in Appendix iBl 

Figure [7] illustrates Type migration in a disk with a power law surface density distribution 
(S oc a~") that decays exponentially on a timescale of 10^ yr. Two of the disks have n = 1; two 
others have n = 1.5. Each pair has an inner edge at 0.1 AU and an outer edge at either 50 AU 
or 100 AU. In all cases, the initial disk mass is 20% of the mass of the 1 Mq central star. The 
migration scales approximately with the amount of disk mass inside a planet's orbit. Therefore 
the effect is strongest for a steep slope in surface density and smaller outer disk radius. While this 
type of migration may not be comparatively strong in magnitude, it may help to alter the orbital 
dynamics in a disk with closely-packed giant planets. 

Type migration pushes planets outward if the disk photoevaporates. It may push planets 
inwards if mass is flowing in toward the star by viscous processes. 



3.1. Migration in a planetesimal disk 



As protoplanetary disks evolve, planetesimals become important to the large-scale orbital 
dynamics of the growing planets. Repeated weak interactions with planetesimals can push a planet 
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radially inwards or outwards. iLin &: Papaloizou I (|l979l ) and iGoldreich fc Tremaine I (|l98d ) first 
quantified aspects of this effect analytically; sinc e then, others have explored it numerically (e.g., 
Malhotralll993 : Ida et al. 2000 : Kirsh et al. 2009 ). Here, we show that our code can perform similar 
calculations. 



Our first illustration of planet migration follows Kirsh et al. (j2009l ). A 2 planet is embedded 
in a disk of planetesimals, each having a mass 1/600*^ of the planet's mass. The disk extends from 
14.5 AU to 35.5 AU and has a surface density S = 1.2 (a/ao)~^ g cm~^ where oq = 25 AU. Initially, 
the planet has e = 0; the planetesimals have e^ms = 

0.001 and 

irms — 0.5 firms 1 where rms' is root- 
mean-squared. The planet stirs nearby planetesimals to large e and ejects some planetesimals from 
the system (Figure [8|). As a reaction to the stirring by the planetesimals, the planet migrates from 
25 AU at 10^ yr to 22 AU at 6 x 10"^ yr, leaving excited planetesimals in its wake. Because the 
planetesimals have finite masses and encounter the planet sporadically, inward migr ation is not 



smoo th (Figure [9]). Our results agree with previous simulations (see Figs. 3 and 4 of iKirsh et al 
2OO9I ). 



Our second demonstration of migration follows the work of iHahn &: Malhotral (| 19991 ). Their 
disk has S oc a^^ and total mass of 100 between 10 AU and 50 AU. The disk is composed of 
1000 planetesimals with Crms = "^irms = 0.01. With the Solar System's four major planets initially 
packed between 5 AU and 23 AU, the system evolves dramatically in time. Dynamical interactions 
between the planets and planetesimals spread the orbits significantly. Jupiter migrates inwards a 
small amount; the outer planets migrate outwards in semimajor axis. After ~ 30 Myr, Neptune 
reaches a ~ 40 A U. Figure 1101 show s the specifics of the migration with our code. The results 
compare well with lHahn Malhotral (|l999l . lower left panel of their Figure 4). 



4. FORMATION OF GAS GIANT PLANETS 

To illustrate how gas giants form in our approach, we adopt the analytic disk model of 
Chambers I (|2009l ) with R^fi = 30 AU, M^fl = 0.01, 0.02, 0.04, or 0.1 Mq, and a = 10"^, lO^^, 10-^ 
or 10~^. We follow the evolution of disk properties (S, M{R), T^, etc) on a grid extending from 
0.2 AU to 10"^ AU. For the coagulation code, we divide the region between 3 AU and 30 AU into 
96 concentric annuli, adopt a dust-to-gas ratio of 0.01, and seed these annuli with planetesimals. 
To set the gas parameters in the coagulation grid, we interpolate physical v ariables in the grid for 



the di sk evolution code. For the bulk properties of planetesimals, we follow iLeinhardt Sz Stewart 
(|2009l ) and adopt parameters for weaker objects with fw = iQh = 2 x 10^ erg g~^ cm''-^, /3fe = —0.4, 



Qg = 0.33 erg g cm /3g = 1.3} (see lKenvon &: Bromlev II2OIOI . and references therein). 



tain (e.g., 


Rice et al. 


2006 


; 


Garaud 


2007: 


Kretke & Lin 


2007 




Brauer et al. 


200 


8; 


Cuzzi et al. 


2008: Laibe et al. 20081: Birnstiel et al. 20091: Chiamr fe Youdin 


2010: Youdin 


2OI0I). To ex- 


plore how the formation of large planetesimals bv the streaming instabilitv fsee Youdin 201C . and 



references therein) might produce gas giants, we consider two extreme possibilities. 



1. Each annulus contains a single large 'seed' planetesimal, r ~ 1000 km, and a large population 
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of 0.1-10 cm pebbles. The gas entrains small planetesimals. The growth rate of the seeds then 
depends on the scale height of the gas, which is set by the viscosity parameter a. Because 
disks with smaller a have smaller vertical scale heights, planets should form more efficiently 
in disks with smaller a. 

2. Each annulus contains only large planetesimals, r ~ 1000 km, with initial e = 10"'*. Aside 

from migration, the gas does not affect larg e planetesimals. However, collision , rates scale in- 

verse ly with the radius of planetesimals (e.g.. lGoldreich. Lithwick. &: Sari II2OO4I : iKenvon &: Bromlev 



20081 ). Thus, planets should form less efficiently than in calculations with a few large seeds 



and a swarm of small pebbles. 



For each of these possibilities, the coagulation code treats the evolution of planetesimals with 
radii of 1 mm to 3600 km. At the start of each calcula tion, we assign a seed to the random number 



gener ator used to assign integral collision rates (e.g., IWetherill k. Stewart 



1993 



Kenvon k. Luu 



19981 ). In this way, otherwise identical starting conditions can lead to very different collision his- 
tories. As long as the gas density is sufficient to entrain particles with r > 1 mm, we distribute 
collision fragments with r < 1 mm into mass bins with r = 1-2 mm. Otherwise, these particles 
are lost to the grid. When individual coagulation particles reach masses of 

1026 g 

(r 3600 km), the code promotes these objects into the A^-body code. After the coagulation code 
initializes a, e, i, and the longitude of periastron, the A'-body code follows the changing masses 
and orbital parameters of promoted objects. Although the A^-body code can treat the evolution 
of A-bodies with any orbital semimajor axis, we assume that objects with a > 10^ — 10^ AU are 
ejected from the system. 

In these calculations, we adopt a single set of parameters for gas accretion onto icy cores. 
Throughout the disk, gas accretion begins when the core mass is 10 and the radius of the 
planet is 0.3 Rj. Gas is accreted cold (r/ = 0.15). Thus, planets are small; mergers of gas giants are 
relatively uncommon. Planets accrete gas fairly slowly, with Mq = 2.0 Mj Myr~^, Mq = 0.1 Mj, and 
am = 2/3. Because planets remain small, they hav e typical luminos i ties 1-2 orders of magnitude 



smal l er than gas giants in other simulations (e.g.. iHubickyj et al. 



2005 



Lissauer &: Stevenson 



2007 



D'Angelo et al. II2OI0I V 



To isolate the importance of disk physics on our results, we ignore gas drag, migration, and 
photoevaporation. For the properties of the disks and planetesimals we consider, growth times are 
generally faster than drag times. In our models, growing protoplanets orbit in a sea of relatively high 
eccentricity planetesimals until they begin chaotic growth. Thus, the conditions required for rapid 



migr ation through a sea of planetesimals are never realized ()Kirsh et al 



2011 



200s 



Bromley Kenvon 



)• 



Many analyses demonstrate that migration through the gas and photoevaporation are im- 
portant pro cesses in arriving at the final distribu t ion of masses arid orb i tal parameters for gas 
giants (e.g., 



Ida &: Lin 



Alexander &: Armitage 



2005 



200s 



Alibert et al. l2005l : iPapaloizou et al.l l2007l : iThommes et al.l I2OO8I : 



Levison et al. II2OIOI . and references therein). In Kenyon & Bromley 



(2011a, in preparation), we show how photoevaporation sets the maximum masses for gas giants as 
a function of a and Mq. By analogy with numerical calculations of planets in disks of planetesimals. 



- 18 - 



Bromley k. Kenvon I (1201 ill speculate that growin g protoplanets are packed too closely to undergo 
type I migration (jWard 1119971 : iTanaka et al.ll2002l ) through the gas until they begin chaotic growth. 
Here, we concentrate on understanding how the growth of planets depends on initial disk proper- 
ties. Kenyon & Bromley (2011b, in preparation) will address how the various types of migration 
change predictions for the masses and orbits of gas giant planets. 

To summarize, our algorithms follow the growth and evolution of planets on three separate 
grids centered roughly on the central star. We derive a solution for the evolution of the gaseous 
disk on a 1-D grid extending from 0.2 AU to 10^ AU. As the most time-consuming part of the 
calculation, we follow the evolution of small solid objects on a smaller, axisymmetric 2-D grid 
extending from 3-30 AU. Once the coagulation code promotes objects, the A^-body code follows 
their trajectories on a 3-D grid extending from the central star to ~ 10^ — 10^ AU. Although objects 
can accrete solid material only from annuli at 3-30 AU, giant planets can accrete gas from 0.2 AU 
to 10^ AU. To save computer time for a large set of calculations with identical starting conditions, 
we halted each calculation at an evolution time of ^ 10 Myr. In future papers, we will investigate 
the long-term stability of planetary systems with gas giant planets. 



4.1. Calculations with Ensembles of Plutos 

We begin with a discussion of planet formation in a disk composed of gas and Pluto-mass 
planetesimals. In these calculations, disk evolution scales with a. Disks with a = 10~^ evolve 
rapidly. For Mdfl = 0.1 Mq, the disk mass declines to 0.07 M© in 0.1 Myr, 0.03 Mq in 1 Myr, 
and 0.008 M© in 10 Myr. Disks with a = 10~^ evolve slowly. For all initial masses, the disk mass 
declines by less than 1% in 10 Myr. 

In disks with ensembles of Plutos, solid objects grow very slowly. In a disk with M^^ o = 
0.1 Mq, it takes only ~ 50-100 yr for the first Pluto-Pluto collision at ~ 3 AU. Due to mutual 
viscous stirring, this first collision occurs at a modest velocity, e ~ 10^^. Thus, this collision yields 
an object with a mass of roughly twice the mass of Pluto and some smaller collision fragments. 
Despite this promising start, it takes ~ 0.1 Myr for this object to accrete another 10 Plutos and 
to reach a mass of roughly 10^^ g. Although these large objects accrete the collision fragments 
efficiently, the fragments have a very small fraction, ~ 0.1-1%, of the initial mass. Thus, accretion 
proceeds solely by the very slow collisions of large objects. 

The dynamical evolution of the growing Plutos is also very slow. At ~ 0.5 Myr, the first 
objects reach masses of ~ 3 x 10^^ g and are promoted into the A^-body code. With typical orbital 
separations of > 0.1 AU, these objects are spaced at intervals of > 10 Hill radii. Thus, they do not 
interact strongly. After ~ 2-3 Myr, the largest objects have masses of 0.1-0.5 Mq and typical orbital 
separations of 0.1-0.2 AU. A few of these are close enough to interact gravitationally. However, the 
interactions are slow. With masses much less than an Earth mass, these objects sometimes reach 
masses of 0.5-1 Mq after 10-20 Myr. These objects will never accrete gas and become gas giant 
planets. 

Figure [TT] shows the mass distributions for these calculations at 3 Myr. Independent of a, the 
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maximum planet mass scales with the initial disk mass as 

oTm©) ■ ^^^^ 

Because the most massive planets result from several rare collisions, the median mass is much 
smaller, ~ 0.05 for disks with M^^q ~ 0.02-0.1 Mq. Disks with smaller initial masses cannot 
produce planets more massive than ~ 0.01 M^. 

For massive disks, the mass distribution is roughly a power law, with a cumulative probability 
p oc m~^. Our results suggest that more massive disks have shallower power laws, with /3 ~ 2 for 
Md,o = 0.1 and /? f» 4 for Mdfl = 0.04 Mq. 

These calculations demonstrate that disks composed of Pluto-mass objects can never become 
gas giant planets. This result is not surprising. For a di sk with surface density S, the accretion time 
is rou ghly t oc RP/T,, where P is the orbital period (e.g. lLissauer 1 ll987l :l iGoldreich. Lithwick. Sari 



2004 '). Thus, large planetesimals grow more slowly than small planetesimals (see also lKenvon &: Bromlev 



20081'). In the limit of large obj ects growing without gravitational focusing, equation (56) of 



Goldreich. Lithwick. &: Sari I (j2004l ) suggests that an ensemble of Plutos can produce 10 objects 



in 200-400 Myr at 3 AU. Our numerical simulations confirm this long growth time. 



4.2. Calculations vi^ith Pebbles and a few Plutos 

Disks composed of a few Plutos and many pebbles evolve rapidly. For a model with M^fi = 
0.1 M0, Rdfl = 30 AU, and a = 10~^, it takes ~ 1.5-1.8 Myr to promote the first object into the 
A^-body code. Within 0.5 Myr, another 10-15 objects reach the promotion mass of 3 x 10^^ g. These 
objects have small atmospheres; they rapidly accrete the remaining pebbles in their feeding zones. 
It typically takes ~ 0.2 Myr for protoplanets to reach masses of 1-2 M^ and another 0.2-0.5 Myr to 
begin to accrete gas from the disk. Within a total evolution time of 2.5-3 Myr, some protoplanets 
grow into gas giant planets. 

The disk viscosity parameter a sets the timescale for the early portions of this evolution. In 
calculations with a = 10^^, the timescale to produce the first A^-body is ~ 2 Myr, much longer 
than the 0.2-0.4 Myr timescale for calculations with a = 10~^ — 10~^. In these models, a sets 
the scale height of the pebbles (equation (j2ip l. Gravitational focusing factors grow with smaller 
scale heights; thus, the Pluto- mass seeds accrete pebbles more rapidly when a is small. Our results 
suggest that planets grow much more rapidly in disks with a < 10~^. 

Stochastic processes are another feature of the evolution. In c alculations with a = 10"*^, 



it takes 2 — 5 x 10^ yr to produce an ensemble of massive oligarchs (IKokubo &: Idalll998l ) with 
Mp > 3 X 10^6 g (Figure [121). For a large set of calculations with identical starting conditions, 
the random timing of oligarch formation and variations in accretion rates for each oligarch lead to 
a broad range in No^max, the maximum number of oligarchs. Often a rapidly accreting oligarch 
prevents a neighboring oligarch from accreting pebbles and reaching the promotion mass. Thus, 
calculations with small No^max tend to have more massive oligarchs than calculations with large 
N 

o,max- 
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Mer gers of protoplanets ar e the final important feature of the growth of gas giant p lanets 



(see also IXhommes et al. I I2OO2I : iKenvon k Bromlev 1 12009 : llda k LinI l20ld : iLi et al. I hom ) . As 



protoplanets grow to masses of 5-10 M^, their orbital separations are typically 10-20 Hill radii. 
Thus, mergers are rare. Once protoplanets start to accrete gas, their separations rapidly become 
smaller than 5-10 Hill radii. Strong dynamical interactions among pairs of protoplanets begin. 
Some interactions scatter protoplanets into t he inner disk at 1-3 AU; others scatter protoplanets 



into the outer disk at > 30-100 AU (see also lMarzari k Weidenschilling II2OO2I : IVeras et al. 



200s 



Raymond et al. I l2009al . 2009b, 2010). Throughout this chaotic growth phase, mergers rapidly 



deplete the number of growing protoplanets (Figure I12|) . Once the number of protoplanets reaches 
~ 5-10, planets have typical separations of > 10 Hill radii. Chaotic growth ends. 

The timescale for chaotic growth varies from one calculation to the next (Figure [T2]) . Some- 
times, protoplanets are widely separated, grow slowly, and reach chaotic growth at late times 
(indigo and blue curves in Figure [T2|) . When several closely packed protoplanets begin to accrete 
gas at roughly the same time, their growth initiates an early phase of chaotic growth (magenta and 
turquoise curves in Figure [T2]) . In nearly all cases, chaotic growth ceases on timescales of 3-10 Myr. 



With growth rates sensitive to a, the mass distributions of planets also depend on a (Figure 
\T3\1 . In calculations with M^fi = 0.1 Mq and R^fi = 30 AU, disks with small a produce a variety 
of gas giants, with masses ranging from 10 to 10-20 Mj. Disks with larger and larger a yield 
smaller and smaller gas giants. For all a, the probability of a given mass is roughly a power law, 
p oc Mp^, with n = 0.2 5-1. Thus, these calculations produce mo re Earth mass planets than Jupiter 
mass planets (see also llda k Lin I l2005l : iMordasini et al. II2009I ). The exponent in the power law 
depends on a: models with smaller a produce broader probability distributions with smaller n. 



For all a, gas giants have a broad range of core masseqj. Although gas accretion begins 
when the core mass is 10 M^, planets continue to accrete solid pebbles. During chaotic growth, 
protoplanets wander through a large range in orbital semimajor axis o, sweeping up (or scattering) 
pebbles and smaller oligarchs along their orbits. Mergers of two gas giants also augment the core 
mass. For disks with M^fl = 0.1 Mq, core masses range from 10 M^ to 100 M^. 

Our calculations suggest that Jupiter mass or larger planets form throughout the disk (Figure 
[T^ . For a < 10~^, Jupiters achieve fairly stable orbits at semimajor axes aj ~ 2-30 AU. Although 
most 10 M0 cores form at a ~ 4 -10 AU, dynamical int eractions scatter protoplanets to o ~ 1-3 AU 



and to a > 15-20 AU (see also iTsiganis et al. I l2005l ) . After the scattering events, these planets 
continue to accrete material and can grow to gas giant planets in their new orbits. 

Planets with smaller masses have broader distributions of semimajor axes than massive gas 
giants. For convenience, we divide lower mass planets into 'Saturns' with masses of 15 M^ to 1 Mj 
and 'super-Earths' with masses of 1-15 M®. For large a ~ 10"^ - 10"^, Saturn -mass planets are 
often the most massive planet in the system and have orbits with as ~ 3-30 AU (Figure [T5]) . When 
a is small, f» 10~^ — 10~^, Saturns are scattered into orbits with as > 100 AU, well outside those 
of Jupiter-mass planets. A few Saturn-mass planets are ejected from the system. Super-Earths 
have even broader distributions of semimajor axis. All super-Earths avoid regions close to gas 



^For simplicity, we assume that all of the solid material in a planet is in the core (see Helled et al. 2008l l 
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giant planets (Figure [T6|) . Although a few super-Earths are scattered into the inner disk, most are 
scattered into the outer disk. Many super-Earths are ejected. 

Lower mass disks produce correspondingly lower mass planets. Figure [T7] shows predicted 
mass distributions for planets in a disk with M^fi = 0.04 Mq. Disks with a > 10^^ fail to produce 
Jupiter mass planets; disks with smaller a rarely produce 3-20 Mj planets. The semimajor axes of 
these planets are much more concentrated towards the central star. The rare Jupiter mass planets 
lie within 10 AU; Saturn mass planets occupy 3-30 AU. Super-Earths occupy the largest range in 
a, with asE ~ 1-100 AU. 

Our results suggest that lower mass disks fail to produce gas giants for any a. In disks with 
Md,o = 0.02 M0 (Figure [T8|l . disks with small a sometimes produce ice giants with Mp 15-30 M®. 
However, most planets in these simulations are super-Earths with Mp ~ 1-15 Mq. Nearly all of 
these planets lie at small semimajor axes, asE ^ 3-10 AU. 

Finally, these calculations produce a diverse set of multi-planet systems (Table [1]) . For any 
Mdfi, disks with smaller a yield a larger frequency of systems with at least two planets. Defining 
Mp^max as the maximum planet mass, planetary systems with smaller Mp^rnax tend to have more 
planets. Systems with 4 Jupiter mass planets are rare: 5-10% among 0.1 Mq disks and less than 1% 
among lower mass disks. Systems with 4-10 super-Earths are relatively common, with frequencies 
of 70% or more among 0.01-0.02 Mq disks. 

5. SUMMARY 

The current version of our hybrid code now includes all of the necessary physics to calculate 
the formation of gas giant planets from an input ensemble of planetesimals in an evolving gaseous 
disk. The new features of the code include 

1. the time evolution of a viscous disk with mass loss from photoevaporation and gas giant 
planet formation, 

2. the atmospheric structure of planets with R ~ 0.001-10 Rq with an algorithm to calculate 
the enhanced radius of a planet and the accretion of small dust grains, 

3. gas accretion onto Earth-mass icy cores, 

4. the time evolution of the radius and luminosity of gas giants, 

5. the gravitational potential of a massive gaseous disk and its impact on the orbits of planets, 
and 

6. the gravitational perturbation of (i) A'"-bodies on a disk of planetesimals and (ii) a disk of 
massive planetesimals on A^-bodies. 



The new hybrid code accurately treats the migration of planets t hrough a disk of massive 
planetesimals. Simulations with isolated planets reproduce the results of lKirsh et al.l (|2009|), where 
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the planet migrates inwa rds through undisturbed planetesimals. Simulations with multiple planets 
reproduce the results of iHahn &: Malhotral (|l999l ). These two examples illustrate the importance 
of the dynamical state of the planetesimals. In the simulations for Fig. [TU[ the initial environment 
of the planetesimals orbiting near Neptune is similar to the initial conditions near the single planet 
in the calculations for Fig. [SHHl Yet, the migration rates in the two cases differ in magnitude and 
sign. Bromley & Kenyon (2011) investigate this phenomenon in more detail. 

To begin to understand the orbital migration of ensembles of gas giants in an evolving, gaseous 
disk, we consider a new 'type 0' migration mechanism. Here, the semimajor axes of planets evolve 
as the disk mass changes. Although the degree of migration is limited by the ratio of the mass of 
the disk to the mass of the host star, modest inward or outward migration is possible. Because 
angular momentum is invariant, the sign of the migration depends on the mode of mass and angular 
momentum loss. In systems that lose mass by disk accretion, migration is inward. Mass loss by 
photoevaporation may lead to outward migr ation. Both modes in ay be important in establishing 
pairs of planets in resonant orbits (e.g., as in lHolman et al. 112010 ). 



To explore the ability of our code to simulate the formation of gas giant planets, we consider 
the evolution of planetesimals in disks with a variety of initial masses and viscosity parameters. 
Disks of Pluto-mass planetesimals cannot form gas giants. For disks with initial masses of 0.1 Mq 
and all values of a, the maximum planet mass is roughly 0.5 Mq. Lower mass disks produce 
substantially lower mass planets. 

Disks composed of 0.1-1 cm pebbles and a few Pluto-mass seeds can produce Jupiter mass 
planets on short timescales. Our results demonstrate that the properties of planetary systems 
depend on the initial disk mass q ^^'i the viscosity parameter a (Figure [T9|) . 



1. More massive disks produce more massive planets on more rapid timescales. In disks with M^fi 
= 0.1 Mq, Jupiter-mass planets are common at 1~2 Myr. In lower mass disks, super-Earths 
and Saturns form in 2-3 Myr. The timescales are similar to the observ ed lifetimes, ~ 1-3 Myr, 
of circumstellar disks around low mass pre-main sequence stars (e.g.. iHaisch. Lada. fc: Lada 



200 ll : iKennedv k Kenvon II2009I : iMamaiek II2009I ) 



^, Jupiters form 
10~^, planetary 



Disks with larger a form lower mass planets. In disks with a = 10~^ — 10 
rarely; Saturns and Super-Earths are common. In disks with a = 10~^ - 
systems with 2 or more Jupiters are common. This result suggests that Jupiters probably 
form more often in disks with 'dead z ones,' regions where the disk viscosity parameter is much 



lower than the rest of the disk (e.g., iMatsumura et al. II2009I ). 



3. The derived distributions of planet mass spans the range of known exoplanets. For planets 
with Mp ^ 0.03-10 Mj, disks with Ma^o = 0.1 M© and a = 10"^ - 10~^ yield a roughly 
power law probability distribution dp/d log Mp oc Mp" with n = 0.20-0.25. Disks with a = 
have a much steeper relation, with n ~ 0.8-0.9. Known exoplanets have an intermediate 
frequency distribution, with n ~ 0.48 (iHoward et al. Il2010l ). Thus, simulations with small 



(large) a produce relatively fewer (more) super-Earths and Saturns than observed in real 
systems. If real protostellar disks have a range of a that includes our small and large a 
regimes, some mix of a can probably 'match' the observations. However, including other 
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physical processes - such as migration and photoevaporatiqn - will change the slop es of the 
mass-frequency relation (see, for example. Ilda &: Lin II2005I : iMordasini et al. II2009I ). Thus, 
detailed comparison with observations is premature. 

4. The derived distribution of semimajor axes is roughly linear in a over 3-30 AU. Scattering 
leads to a few gas giants at 1-3 AU and at 30-100 AU. Without a treatment of migration, our 
predicted distribution of semimajor axes cannot hope to match observations. However, our 
ability to scatter gas giants out to 30-100 AU m ay allow core accretion models to produce 



systems similar to HR 8799 (jMarois et al. II2008I ) 



5. These calculations often produce planetary systems with 2-10 planets. Planetary systems 
with multiple planets usually contain super-Earths, sometimes contain Saturns, and rarely 
contain Jupiters. Current observations suggest that r oughly 10% of all planetary systems 



have two or more gas giant and/or ice giant planets (jGould et al. Il2010l . also statistics at 



exoplanet.eu and exoplanets.org). However, some planetary systems may have as many as 



7 super-Earth or ice giant planets (e.g., HD10180; iLovis et al. Il20ld ). Given current poor 



constraints on initial disk mass and viscosity in protostellar disks, our results can 'match' 
both of these observations. 

6. Jupiter mass planets in these calculations have core masses of 10-100 Mq. Once icy 10 
cores start to accrete gas, they continue to accrete solids from the disk. Mergers with other 
gas giants also augment the core mass. When mergers between growing protoplanets are rare, 
planets have modest core masses, 10-20 M^. When mergers between icy cores are common, 
core masses approach 100 M®. Multiple mergers of growing protoplane ts may account for the 
large heavy element abundances in HD 149026b and other exoplanets (jBurrows et al. 1120071 : 



Fortnev et al. Il2009l l 



These results are encouraging. Our simple set of initial calculations accounts for the observed 
range of planet masses and for some of the observed range of semimajor axes. The derived power law 
slopes for the frequency of planet masses bracket the observations. Compared with real planetary 
systems, the calculations probably produce too many multi-planet systems. 

Adding more realism (e.g., migration, photoevaporation, etc) to our simulations will allow us to 
improve our understanding of the processes that lead to ice/gas giant planet formation. Integrating 
orbits past ~ 10 Myr will enable robust comparisons between the properties of real and simulated 
planetary systems. Together with improved observations of protoplanetary disks and better data 
for exoplanets, numerical simulations like ours should lead to clear tests for the core accretion 
theory of planet formation. 
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A. The gravitational potential of an unperturbed disk 

If the gaseous protostellar disk is extended, thin, and axisymmetric, our code can treat the 
gravitational effect of the disk on the planets. We estimate the gravitational potential produced by 
a geometrically thin, axisymmetric disk in two ways. The first method is based on an analytical 
inversion of the Poisson equation from an expansion of the surface density S into Bessel functions. 
The results are simple expressions for the potential, valid for 2D disks of large spatial extent and a 
surface density that is a pure power-law in orbital distance. The second method employs brute force, 
using a direct numerical integration of the mass density over a geometrically thin, axisymmetric 3D 
disk. This approach gives the potential associated with an arbitrary surface density profile, with 
the advantage of generalizability at the expense of computational load. 

The analytical approach for estimating the disk potential as it acts on a unit-mass particle 
takes advantage of solut ions to the Laplace equatio n that are linear combinations of Bessel functions 



of the first kind, Jo{x) (jBinnev fc Tremaine 



19871 . §2.6) 



POO 

^diR,z)= S{k)Jo{kR)e-''\'kk, (Al) 
Jo 

where S{k) is a density function. The exponential with the discontinuous first derivative at z = 
generates a 6 function in the 3D density at the disk midplane upon application of the Laplace oper- 
ator. Gauss's law relates this potential to the surface density. A closed, pillbox-shaped surface with 
a unit-area cross-section that straddles the disk midplane has a gravitational field fiux —2d^d/dz, 
and encloses a mass S. Gauss's law sets the fiux equal to AttG times the enclosed mass, hence 

= / S{k)MkR)kdk. (A2) 

2vrG Jo 

The integral is a Hankel transform (to within a multiplicative constant), and we identify S{k) and 
E(i?) as transform pairs, with 

S{k) = -2itG / i:{R)Jo{kR)RdR. (A3) 
Jo 

We now assume that the surface density S(i?) is a power law in R with an arbitrarily large 
radial extent, i.e., -Rin — )• and Rout — > oo. In this case, 

S(a) = So(i?o/i?)-" (A4) 

and 

poo 

Sik) = -2TTGEoR^k''-^ / q^-''Jo{q)dq . (A5) 

Jo 

We impose the condition that l/2<n<2so that the mass in the disk is finite inside a finite 
radius, with the lower limit set so that the integral in equation IA5I is well-defined. This integral is 
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of the form 



10 



Inili) = I x^'Jn{x)dx (A6) 



= 2^ ^ J > (A7) 



1 



(A8) 



S{k) « -27rGSoiiSfc''"^-fo(2 - n) . (A9) 
Similarly, the potential in the disk midplane is 

= -27rGSoE (^-|-) /o(l-n)/o(n-2), (AlO) 

although this expression is formally valid only for 1 < n < 2. The radial component of the 
acceleration, a^, is not subject to this formal restriction on n; it obtains from the derivative of 
with respect to R, along with the identity dJo{x)/dx = —Ji{x): 

ar = -27rGSo/o(l - n)/i(l - n) f — ^ , (All) 

which is exact in the midplane of an infinitely extended disk. 

To illustrate results from the above analysis, we give the following examples of the potential 
and radial acceleration for a few power-law surface density profiles: 



0.6 : = 2.87rGEoi?[]-6i?°-^ = -l.l7rGEo(i?/-Ro) 



-0.6 



n = 1.0 ^d = 27rGSoi?o log(-R) + const = -2txGT.qRq/R (A12) 

n = 1.5 = -8.87rGSoi?(^-^i?-°-^ = -4.4^GSo(i?/-Ro)"^-^ 

As long as the orbital distance R approximately satisfies iijn < O.li? and lOi? < i^outi these 
expressions hold for disks of finite extent. 

As an alternative to the analytical method, which is limited in terms of the functional form of 
we take a numerical approach to directly integrate the mass in the disk to get the gravitational 
potential: 

jRin J-^ \fiTl_ + z^ + h^ + R^ - 2r^Rcos{ip) 

where position r relative to the host star is in cylindrical coordinates, {r±,ip,z). The variable h is 
a softening parameter, which formally eliminates the density singularity in the disk midplane by 
making the 2D disk behave as if it has some finite thickness. From the Poisson equation, we find 
that the effective vertical density profile is approximately Gaussian, with a scale height of a ~ /i. 
Here we take h = 0.01 AU; in what follows, the results do not strongly depend on this choice. 

The next step is to integrate eq. (|A13P numerically; the radial integration can be solved ana- 
lytically, but integrating over the angular coordinate if requires a Chebyshev-Legendre quadrature 



-26- 



scheme. The number of sample points n in our quadrature scheme is variable; we increase n by a fac- 
tor of four successively until the relative change in the result from one iteration to the next is below 
an error tolerance of 10~^. For the integrands encountered here, the result is much more accurate 
than the error tolerance suggests, because the quadrature scheme typically converges quadratically 
or better. 

Compared to with evaluating the 1/r potential from the central star, this method is compu- 
tationally expensive. To alleviate this problem, we limit the calculations to the disk midplane — a 
reasonable approximation since the orbital inclination of planets not entrained in the gas is typ- 
ically less than the disk scale height. Then we evaluate the potential at a set of logarithmically 
spaced points in orbital separation and interpolate with cubic splines. The computational load 
then reduces to 0(10) arithmetic operations. 



B. Time- varying disk mass ("type 0" migration) 

When the disk mass is removed from the planetary system by photoevaporation or some other 
slow process, the time- varying disk mass can produce a radial migration of a planet's orbit. Angular 
momentum is conserved. For example, the product of the semimajor axis a and the central star 
mass M remains constant when the star itself loses mass. For low-eccentricity orbits, the conserved 
quantity is equal to a^a^, where ar is the acceleration on the body from both the star and gaseous 
disk. If the disk mass — but not the stellar mass — varies in time, 

= l + 2.Eo/M/o(l-n)/.(l-n)»(Of-- 

^ ^l + M27rSo/M/o(l-n)7i(l-n)aga(t)2-ne-*/^ ^ ^ 

« [l-h27rEo/M/o(l-n)/i(l-n)aSa(0)2-^] (i > r) . (B2) 

Setting n = 1 and taking the limit of t ^ oc, the maximum change in orbital distance for a planet 
at 1 AU in an evaporating disk with initial surface density SO = 3000 g/cm^ is roughly 1%. If a 
disk with the same mass is more centrally concentrated (n = 1.5), this number is roughly an order 
of magnitude higher. For planets farther out in the disk, the maximum change is much larger, ~ 
5% at 10 AU and ~ 10% at 30 AU. 
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Table 1. Frequency of Multi-Planet Systems^ 
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^Fraction of planets in systems where the most massive planet is a Jupiter ('J' 
in Mass Key column), a Saturn ('S'), or a super-Earth ('SE'). For example, in the 
first row, disks with M^^ = 0.1 M© and a = 10~^ produce Saturn or lower mass 
planets; the frequency of systems with N Saturn- mass planets is 30% [N = 0), 
16% (iV = 1), 30% (iV = 2), 16% (iV = 3), 5% {N = 4) and 0% (iV > 5). 
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Fig. 1. — Time evolution of the surface density of a gaseous disk (M^ q = 0.04 M0, Rii q = 10 AU, 
a = 10~^) surrou nding a 1 Mq star. Dashed hues show results for the analytic disk model of 



Chambers I (|2009l ): solid lines show results for our numerical solution of the diffusion equation. 
Despite small differences in the initial conditions, the numerical solution tracks the analytic model. 




Fig. 2. — Time evolution of the disk mass (upper panel) and disk accretion rate onto the central 
star (lower panel) for the analytic and numerical solutions in Figure 1. 
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Fig. 3. — As in Figure [2] for disks with photoevaporative winds. The legend in the upper panel 
indicates the mass loss rate of the wind in units of yr~^. 
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1 ' 1 ' 1 ' 1 r 




log Time (yr) 

Fig. 4. — Luminosity evolution for a gas giant planet using the prescription in equations (f26|) - ([28]) . 
Starting with an initial mass of 10 and an initial radius of 1 Rj, the planet accretes gas at a 
rate of 10~^ Mj yr~^ until it reaches a mass of 1 Mj. Solid curves show the evolution for various 
rj. When r] is small; the planet accretes cold material and has a small peak luminosity. When rj 
is large, the planet accretes hotter material, expands, and has a large peak luminosity. At late 
times, all curves converge on a single evolutionary track which yields a Jupiter-radius planet at 4-6 
Gyr. The dashed curve illustrates how the evolution changes for planets with initial radii of 0.33 
Rj (comparable to the radius of Neptune). 
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Fig. 5. — As in Figure H] for the radius of the planet. Planets that accrete hotter gas grow to larger 
radii. 
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Fig. 6. — The evolution of (e^ + i'^Y ^'^ (in units of Hill radii) of a particle disk with two embedded 
planets, as in iKokubo Ida I (|l995l ). The abscissa gives the relative orbital separation in terms 
of the difference between the semimajor axis a and a reference radius uq = 1 AU. The histograms 
correspond to the disk particles; the symbols represent the two planets. The colors and symbol styles 
distinguish simulations: black histograms and circles indicate fully interacting massive particles, 
blue lines and triangles show data from a massive disk of swarm particles with no self-gravity, while 
magenta histograms and squares correspond to a massless disk composed of tracers. 
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Fig. 7. — Scmimajor axis as a function of time in Type migration. Planets move in response to 
a decrease in disk mass, for four different massive disks, with surface density S ~ a^'\ and outer 
disk radius aout as labeled. The initial mass of each disk is 20% of the stellar mass (1 Mq); the disk 
mass decays exponentially on a timescale of 10*^ yr. The greatest migration is for the n = 1.5 disk 
with an inner edge of 50 AU (gray dashed line), which has the most mass initially concentrated 
within the orbits of the planets. 
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Fig. 8. — The evolution of eccentricity in units of Hill radii as a function of semimajor axis during 
planet migration. The three panels show snapshots at the in dicated times as a 2 planet travels 
through a sea of planetesimals, similar to lKirsh et al.l ( 2003 ) — S ~ 1/a, normalized to 1.2 g cm^^ 
at 25 AU, and extending from 14.5 AU to 35.5 AU. The blob of particles in each panel (~ 25 AU 
at 10^ yr, 24 AU at 3 x 10^ yr, and 22 AU at 6 x 10^ yr) indicates the location of the planet. 
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Fig. 9. — The evolution of the semimajor axis of a planet in the planetesimal disk from Figure El 
The initially flat migration profile stems from "inertia" associated with planetesimals in the co- 
orbital zone (the blobs of particles in Fig. [51 which diffuse away in time) . After the planet scatters 
these particles out of its orbit, fast migration begins. 
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Fig. 10. — Migration of the major planets in a planetesimal disk. The planets, with their present- 
day masses, evolve from an initially compact configuration as a result of interactions with 1000 
planetesimals, each with a mass of 0.1 M®, uniformly distributed in semimajor axis from 10 AU to 
50 AU. In contrast to the behavior of the single planet in Fig. [9l interactions bet ween planetesimals 
and m ultiple planets cause the outward migration of Neptune and Uranus (see iHahn &: Malhotra 
1993) • 
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Fig. 11. — Mass distributions at 3 Myr for calculations with ensembles of Pluto-mass planetesimals. 
The legend indicates the initial disk mass in Mq. The maximum planet mass scales with the disk 
mass. For planets with masses exceeding 10 M^, the cumulative probability is roughly a power 
law, p oc Mp^ . Less massive disks have steeper probability distributions. Disks with initial masses 
Md.o < 0.01 Mq do not produce 10 M® planets. 
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log Time (yr) 

Fig. 12. — Time evolution of Nqu, the number of oligarchs with Mp > 3 x 10^^, in models with 
Mdfi = 0.1 Mq, RfiQ = 30 AU, and a = 10^^. In these simulations, N^h rises as Pluto mass seeds 
rapidly accrete tiny pebbles, reaches a plateau when these seeds begin to accrete gas, and then falls 
as gas giant planets merge with other gas giants and smaller oligarchs. As gas giants achieve stable 
orbits. Noli stabilizes. 
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Fig. 13. — Mass distributions at 3 Myr for disks with Mdfi = 0.1 Mq, Rdfi = 30 AU, and a as 
indicated in the legend. Disks with smaher a produce more massive gas giants. The median mass 
of gas giants ranges from ~ 20 Mq for a = 10~^ to ~ 0.25 Mj for a = 10"^. 
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Fig. 14. — As in Fig. [13] for the semimajor axes of planets with masses larger than 1 Mj. Disks 
with smaller a produce gas giants over a broader range of semimajor axes. The median semimajor 
axis ranges from ~ 6 AU for a = 10^'^ to ~ 8 AU for a = 10~^. Disks with a = 10^^ do not make 
massive gas giant planets. 
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Fig. 15. — As in Fig. [TUfor the semimajor axes of planets with masses between 15 and 1 Mj. 
Disks with smaher a produce Neptune to Jupiter mass planets at larger semimajor axes. At 3 Myr, 
the median semimajor axis ranges from ~ 6 AU for a = 10~^ to ~ 50 AU for a = 10~^. 
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Fig. 16. — As in Fig. [T3] for the semimajor axes of planets with masses 1-15 Mq. Disks with 
smaller a produce super-Earth mass planets over a broader range of semimajor axes. Super-Earths 
avoid regions near Jupiter mass planets. At 3 Myr, the median semimajor axis ranges from ~ 5 AU 
for a = 10-2 to « 60 AU for a = 10"^ 
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Fig. 17. — Mass distributions at 3 Myr for disks with Mdfi = 0.04 Mq, Rdfl = 30 AU, and a as 
indicated in the legend. Disks with smaller a produce more massive gas giants. The median mass 
of planets ranges from ~ 1.5 for a = to ~ 20 for a = 10~^. 




Fig. 18. — As in Figure [T71 for disks with M^fi = 0.02 Mq. Note the change of scale for the abscissa 
compared to other figures. The median mass ranges from ^ 1 Me for a = 10"^ to ^ 6 for 
a = 10-5. 
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Fig. 19. — Outcomes of planet formation simulations. The regions indicate the maximum masses 
of planets as a function of the initial disk mass, M^fi, and the disk viscosity parameter, a. 



